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ABSTRACT 

To better model the efficient production of cosmic rays (CRs) in supernova remnants (SNRs) with the associ- 
ated coupling between C R production and SNR dynamics, we have generalized an existing cr-hydro-NEI code 
(i.e.. lEllison et al.|[2012l) to include the following processes: (1) an explicit calculation of the upstream precur- 
sor structure including the position dependent flow speed, density, temperature, and magnetic field strength; 
(2) a momentum and space dependent CR diffusion coefficient; (3) an explicit calculation of magnetic field 
amplification (MFA); (4) calculation of the maximum CR momentum using the amplified magnetic field; (5) 
a finite Alfven speed for the particle scattering centers; and (6) the ability to accelerate a superthermal seed 
population of CRs as well as the ambient thermal plasma. While a great deal of work has been done modeling 
SNRs, most work has concentrated on either the continuum emission from relativistic electrons or ions, or the 
thermal emission from the shock heated plasma. Our generalized code combines these elements and describes 
the interplay between CR production and SNR evolution, including the nonlinear coupling of efficient diffu- 
sive shock acceleration (DSA), based mainly on the work of P. Blasi and co-workers, and a non-equilibrium 
ionization (NEI) calculation of thermal X-ray line emission. We believe our generalized model will provide 
a consistent modeling platform for SNRs, including those interacting with molecular clouds, and improve the 
interpretation of current and future observations, including the high-quality spectra expected from Astro-H. 
SNR RX J1713.7-3946 is modeled as an example. 

Subject headings: acceleration of particles, shock waves, ISM: cosmic rays, ISM: supernova remnants, mag- 
netic fields, turbulence 



1. INTRODUCTION 

In this paper we report on an extensive generalization of 
a cr-hydro-NEI code (e.g.. JEUison et al.ll2007l: iPatnaude et alJ 
l2009l l2OT0l Ellison etal.ll20ia 120121) that has been used 
to model cosmic-ray (CR) production in evolving supernova 
remnants (SNRs). This code uses a semi-analytic, nonlin- 
ear diffusive shock acceleration (DSA) model to calculate the 
CR production at the remnant forward shock (FS) and deter- 
mines the broadband continuum radiation from CR electrons 
and ions consistently with the thermal X-ray line emission 
from the shock-heated plasma. The feedback effects from ef- 
ficient CR production influence the SNR dynamics and the 
non-equilibrium ionization (NEI) calculation of the thermal 
emission. Until now, th is code has used the nonUnear D SA 
calculation described in lBlasi. Gabic i. & V annonil (I2005h . In 
a series of papers since then (e.g..lCaprioli et al.l2009l and ref- 
erences therein), Blasi and co-workers have presented a num- 
ber of important generalizations to their semi-analytic model 
and we have included many of these, as well as additional 
effects, in an updated version of the cr-hydro-NEI code. 

We note that in addition to the work of Blasi and co- 
workers, a great deal of work on nonlinear DSA has been 
done by a large number of researchers in the last 30 years 
(for reviews see Drurv 1983; Blandford & Ostriker' 119781 
Jones & Ellison 1991; Malkov & Dmry 2001; Schu re et al.l 
20121) . A good descriptio n of some of the nu merical 
techniques used is given in ICaprioli et al.l (1201 Obi) where 
a distinction is made between fully numerical time- 
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depen den t solutions (e ^.. | B elll 119871; iBerezhko & VHlk 
1997[ iKang et all I2009E IZirakashviU & Aharonian 
20 IC), station ary Mo nte Carlo simulations (e.g., 
Jones & EllisonI 119911: IVladimirov. Bvkov. & EUisonl 



20091). and quasi-stationa r y sem i -analytic solutions (e.g., 
Malkov. Diamond. & VoUd l2000t iBlasi. Amato. & CaprioUl 



ft 



2007t ICaprioli. Amato. & Blasil boiOal) . In addition to 
these techniques there are direct particle-in-cell (PIC) 
plasma simulations that, in principle, can model all 
aspects of shock formation, particle injection and accel- 
eratio n, and the self-generation of magnetic turbulence 
(e.g., iHoshino & Shimadal l2002t jRiauelme & SpitkovskyI 



20 id I2OIII; iGargate & Spitkovskvf 12012b (see, however 
Vladimirov. Bykov. & EllisonI 120081 for a discussion of 
the difficulties in applying PIC simulations to SNRs). We 
have chosen to adapt the semi-analytic Blasi model to our 
cr-hydro-NEI code because this model contains most of the 
essential nonlinear effects we need, is computationally fast, 
and is fully documented in published papers. 

The specific generalizations we include are: (i) the ex- 
plicit calculation of the upstream precursor structure includ- 
ing the position dependent flow speed, density, temperature, 
and magnetic field strength; (ii) a momentum and space de- 
pendent CR diffusion coefficient; (iii) an explicit calculation 
of magnetic field amplification (MFA) that replaces our pre- 
vious ad hoc parameterization; (iv) the calculation of the 
maximum CR momentum, Pniax, using the amplified mag- 
netic field; (v) a finite Alfven speed for the particle scattering 
centers; (vi) the ability to accelerate a superthermal seed pop- 
ulation of CRs as well as the ambient thermal plasma; (vii) a 
calculation of the break-momentum and high-energy cutoff of 
the CR distribution due to ion-neutral damping when the FS 
interacts with a partially ionized medium; and (viii) Coulomb 



losses for low energy but still superthermal CRs. 

Many assumptions and approximations are required for 
these generalizations and we detail them here. We compare 
results from our new generalized code to previous work mod- 
eling SNR RX J 17 13.7-3946 (henceforth SNR J 17 13) and 
point out where significant changes have occurred. While 
our generalized model sh ows differences from the descrip- 
tion of SNR J 1 7 1 3 given in lEllison et alj (l2012h . our most im- 
portant conclusion, that the GeV-TeV emission is dominated 
by inverse-Compton (IC) emission from relativistic electrons 
rather than pion-decay emission from relativistic ions, re- 
mains robust. The fit to SNR J 17 13 from our more complete 
and consistent model is accomplishe d with a modest vari ation 
in parameters from those detailed in lEllison et al.l (12012b . 

In addition to fitting SNR J1713, we verify that our adap- 
tation of the Blasi model matches the published work of the 
Blasi group, where appropriate, and show cases where our 
model differs from the implementation presented by Blasi 
and co-workers, detailing why these differences occur The 
most important difference is that we have coupled the up- 
dated steady-state model of Blasi and co-workers to the time- 
dependent, SNR hydrodynamics Q Other differences concern 
the treatment of the CRs escaping at the free escape boundary 
(FEB) and how Pmax is determined. 

An alternative technique using a Monte Carlo simula- 
tion of nonlinear (NL) PSA (e.g iJones & ElUsonl 119911: 
IVladimirov. Bvkov. & Ellisonll2009b makes a different set of 
assumptions from the semi-analytic calculations and we have 
compared our results against this model as well. In this case, 
significant differences occur particularly concerning thermal 
particle injection. There are less pronounced differences in 
the maximum particle momentum, Pmax- 

While our cr-hydro-NEI code includes the influence of par- 
ticle acceleration on the evolving SNR, the calculation of 
the thermal X-ray line emission coupled to the CR produc- 
tion, and the production and propagation of escaping CRs 
to nearby mass concentrations, it remains a spherically sym- 
metric model that doesn't attempt to describe the small-scale 
structure seen in most SNRs. However, despite these approx- 
imations, we believe our generalized model is currently the 
most complete description of a young SNR in terms of the 
broad-band emission from radio to TeV from a segment of the 
remnant large enough to safely average over the small-scale 
structure. Future plans for the cr-hydro-NEI code include cal- 
culating the radiation emitted at the reverse shock and model- 
ing the interaction of CRs with dense molecular clouds. 

2. MODEL 

Supernova remnants are complex phenomena and any 
model with a claim to self-consistency between the rem- 
nant dynamics and the nonlinear production of CRs will 
be complicated. The time-dependent SNR hydrodynam- 
ics must be coupled to the CR production since the ac- 
celeration of CRs changes the ratio of specific heats for 
the shocked plasma and drains energy from the plasma 

' Some justification for using a steady-state DSA calculation in an evolv- 
ing SNR is necessary. Basically, we assume that as long as the dynamic 
time-scale of the SNR is long compared to the acceleration time to Pmax. 
tacc(Pinax), the Steady-State approximation will be valid. In all of the ex- 
amples we show here, except for very early times (i.e., t < 50 yr), we have 
tacc(pmax) <C isNR- We note that the quasi-steady-state model we use 
is also in good agreem ent with time-dependen t models that show "CR dilu- 
tion" (see, for example IBerezhko et al.lll996alibtlEnTson & Bvkovll20Tll and 
references therein). 



as the highest energy CRs escape from the FS. The 
DSA of CRs is difficult to model in a static config- 
uration, particularly when the nonlinear shock structure, 
MFA, and the self-consistent determination of the diffusion 
coefficient is included (e.g., [V ladimir ov. Bykov. & EUisonl 
I2008t iBykov. Osipov. & EUison 201 Ih). An evolvin g SNR 



IE 



adds additional complications (e.g., Berezhko et al.lll996at 
IZirakashvili & Ptuskin||2008l: iBvkov et al.ll201 lal) . 

Given these complications, it is tempting to suggest that 
a particle-in-cell (PIC) simulation, rather than an analytic 
or semi-analytic calculation, is the best approach (e.g., 
iRiquelme & Spitkovskyll2010l) . In principle, a PIC simula- 
tion including all of the coupled phenomena could be done 
but in practice, for the immediate future at least, the dy- 
namic range for particle energies and diffusion lengths pre- 
clude such a direct approach for th e nonrelativistic shocks in 
SNRs (see IVladimirov et a l.' 2008). On the other hand, fuUy 
anal ytic self-similar (e.g.. Chevalier 19 821) or parameter ized 
(e.g.. lTruelove & McKeell 999; Morlino & Capriolil20Ilh ap- 
proaches lack the ability to include nonlinear effects or MFA 
self-consistently. Monte Carlo techniques have been used 
extensively to model nonlinear DSA including MFA (e.g., 
IVladimirov et al. 2008, 2009). However, these techniques are 
computationally intensive and intrinsically steady-state and 
are not well-suited for inclusion in a SNR hydro simulation. 
We beheve the semi-analytic approach developed here pro- 
vides a compromise that is accurate enough to model obser- 
vations of an evolving SNR and extract critical information on 
the supernova and environmental properties until more com- 
plete numerical simulations are "large" enough to fully model 
an evolving SNR. 

When relativistic CRs are produced efficiently by the SNR 
shock, the pressure to energy density ratio in the shocked 
plasma decreases. In addition, some fraction of the highest 
energy CRs will escape from the SNR further reducing the 
pressure driving the FS. These effects result in a shock com- 
pression ratio larger than would be the case with test-particle 
CR production and impact the SNR evolution. The magnetic 
field amplification further modifies the hydrodynamics by in- 
creasing the magnetic pressure downstream from the shock. 
Our model couples the hydrodynamics, DSA, and MFA al- 
lowing a self-consistent calculation of the CR spectra and re- 
sulting continuum radiation, as well as the thermal X-ray line 
emission in the interaction region between the FS and CD. 
The changes produced in the hydrodynamics by NL DSA 
modify the ionization state of the shocked plasma and these 
changes are self-consistently included in our non-equilibrium 
ionization calculation of the thermal X-ray emission. The de- 
tails of our generalized code are given below. 

2.1. Supernova remnant and shock hydrodynamics 

As in previous work (e.g.. lEllison et al.ll2012h . the SNR is 
modeled with a spherically symmetric hy dro simulation based 
on VH-1 (e.g.. iBlondin & Enisonll200l PI The standard hy- 
drodynamic equations are modified when efficient DSA oc- 
curs through a change in the equation of state from the influ- 
ence of relativistic CRs, energy loss from escaping CRs, and 
magnetic pressure. The equation of state used in the hydro 

^ The basic hydro simulation VH- 1 is open to the public and can be found 
at |http : / /w onka .physics .ncsu .edu/pub/VH-1/ This web site 
includes a user guide and full background information on the code. We em- 
phasize, however, that while the hydrodynamics are essentially the same as 
in this basic code, the ve rsion we employ has been modified extensively to 
include NL DSA (see Ellison et al. 201^ and references therein). 



equations is modified by using the ratio of specific heats cal- 
culated from the particle distribution function so that the ratio 
lies between 5/3 and 4/3 depending on the fraction of nonrel- 
ativistic to relativistic particle pressure. The SNR simulation 
provides the shock speed, sonic and Alfven Mach numbers, 
and other quantities necessary for the shock acceleration cal- 
culation which is performed as follows. 

From the conservation of momentum flux, we can write the 
following equation in the shock rest-frame: 



Po^l- 



Pgfl + PCR,Q + Pwfi = 

^2 



pixMxY + Pgix) + PcRix) + PUx) . (1) 



Throughout this paper, we adopt a coordinate system in the 
shock rest-frame such that the subshock is located at x = 0, 
and the subscripts '0', T and '2' refer to regions far up- 
stream (x = — oo), immediately upstream [x = 0^) and 
immediately downstream (x = 0^) from the subshock, re- 
spectively. Equation ^ implicitly assumes that the shock 
is plane and we will always assume that the shock precur- 
sor is a small enough fraction of the shock r adius for this to 
be a valid approximation (see, for example, iBerezhko et al.l 
Il996at iPtuskin et alJl2010l for models that explicitly include 
the spherical symmetry and CR dilution). Equation dU fur- 
ther assumes that the shock obliquity is unimportant, i.e., that 
the upstream magnetic field is nearly parallel to the shock 
normal. Essentially, we make the assumption that the mag- 
netic field is tangled enough so that all of the effects of the 
magnetic field can be accounted for through a scalar mag- 
netic pressure. In equation ([l), po is the far upstream den- 
sity including helium if it is present, 7^ = 5/3 is the ratio 
of specific heats for the thermal gas, A/q is the sonic Mach 
number, Pgo = poUQ/{jgMo) is the ambient gas pressure, 
-PcR.o is the pressure from any pre-existing seed CR proton 
population, Pcb.{x) includes the pressure from the acceler- 
ated protons, Pu, is the pressure from any pre-existing am- 
bient magnetic turbulence, and P^, (x) includes the pressure 
from the self-generated magnetic turbulence through resonant 
CR streaming instability. 

We include precursor heating from th e damping of turbu- 
lent w aves and adopt equation (50) from IBerezhko & EllisonI 
(Il999h . such that the gas pressure at any a: < in the precur- 
sor region is given by 



Pgix) - Pg^oUixy^ X 



M^ 



l + /damp(73-l)T7^(l-C/(a:P) 



(2) 



where U{x) = u{x)/uq and Ma,c« = \/Maj)Ma^ = 

^JuoUl/{vAfiVA,l)^ Here we have assumed that the speed 
of the Alfven waves va is unimportant compared to the shock 
speed uq, which is always the case for the models described 
in this paper A model of wave-damping such that the wave- 
growth rate a{x) is proportional to the damping rate r{x) is 
used, so fliat r(a;) = /dampO-(x), where /damp G [0, 1] is a 
free parameter of the model. 

For the magnetic pre s sure, w e refer to the quasi-linear 
model in ICaprioli et al.l (l2010ah and express the magnetic 

' Note that in many formulations M^ q i** assumed in equation (2) rather 
than A/4 ;,ff but for most cases where the Alfven Mach is large, this is not 
significant. 



pressure in the precursor region in terms of U{x) as 
dBix)"^ 



PM) 



Stt 

AR2 

"-°ISM 



(1 - /damp)/30^to fl- U {x) 



AMa.O 



U{x)y^ 



,(3) 



where 6Bism is the background turbulent field strength in the 
far upstream medium. For simplicity, we do not consider the 
presence of pre-existing magnetic turbulence in the ISM and 
assume SBism = (or P^.o = 0) for all models presented 
in this paper The total amplified i3-field in the precursor is 
hence 



B{x) ^ \6B{xf + Bl 



/87r(l - /damp)yOoMo (^- U{x) 



4Ma,o 



C/(a;)3/2 



B'o 



(4) 



since we assume that the generated turbulence is purely 
Alfvenic and transverse to the background field Bq. For mod- 
els in which MFA is not included, Puj (x) is set to zero every- 
where. 
The CR pressure Pen is given by 



Pcr{x) = PcR.oUix)--'^''^-^ 



+ 



4^ /-Pmax 

-IT / P^v{p)f{x,p)dp , (5) 

■^ "'Pinj 

where f{x,p) is the phase-space distribution of the acceler- 
ated particles and 7cR(a:^) = 1 + Pcr{x) / Ecr{x) is the ra- 
tio of specific heats for the CR particles. The production of 
f{x,p) via nonlinear DSA is described in Section IZ2] Fi- 
nally, we can rearrange equation ([!} in a dimensionless form 
by dividing it with pqUq: 



^CRix) = 1 - U{x) + 



1 



XCR 



+ 



./damp(7g - 1) 



79MS 



7gMA,cff 

/damp f l- U{xY 
[/(x)3/2 



(l-C/(x)-^«) + 



4Ma.o 



(6) 



where CcR(a;) = PcR{x)/{paul) and xcR = ^CR,o/^g,o- 
The third term on the right-hand-side (RHS) describes the 
adiabatic compression and wave-heating in the precursor, the 
fourth term on the RHS describes the pre-shock compression 
of the incoming pre-existing non-thermal particles, and the 
last term is t he pressure froni the se lf-generated turbulence. 

Following ICaprioli et al.l (l20I0ah . the total and subshock 
compression ratios experienced by the background fluid, i?tot 
and Rsnh respectively, are related through the following equa- 
tion: 



R 



.7s + l . 



2z?Ts 



M^R 



0-"^sub 



Ab = 



_ -Ptu.l 



n,l 



(7g + 1) - Rsnhilg - 1) 

(1 + Ab)(1+Ath) 
1 + i?sub ( — - 1 



7g 
M2 



ATH^/damp(7, " 1)t7^(1 ^ U^^ 



(7) 



where As and Ath give the corrections from the magnetic 
jump condition with self-generated turbulence, and from pre- 
cursor heating, respectively. Equipped with the compression 
ratios, one can readily write down the gas temperatures in the 
upstream precursor T{x) and immediately downstream from 
the shock, T2: 



r(x)=To 



Pgi^) 



P. 



3,0 



Uix) 



= ToU{x) 



1-79 



M2 



and 



MA,eS 



Pn 



, (8) 



T2=Ti ( '-^\R-^ 



'sub 



= T, 



n,i 



(7a + l)fls„b-(7s-l) + (7a-l)(fl.ub-l)A 
(7g + l-(7s-l)-R=ub)flsub 



/^^Pt^-^lu _ ^\2 



n.l 



(^sub — 1)^ 



^ 7i,Mo^(l^/d„.p)(l-C/i^)(fl,„b-l)" 

4A/A,o[l + /da„.p(7s-l)(M2/MA,eft)(l-C/7'')]C/i^''~'"' ' 

(9) 

We note that the hydro code now includes the CR precur- 
sor self-consistently producing a smooth transition from the 
far upstream values of density, pressure, temperature, flow 
speed, and B-field to the subshock values. This has been in- 
cluded for consistency but does not significantly modify the 
hydrodynamics because the precursor drops sharply upstream 
from the subshock, as shown in Figure [T]^ However, the pres- 
ence of a precursor is an important prediction and diagnostic 
for NL DSA and our calculation will be important for future 
determinations of the precursor emission]] 

2.2. Diffusive Shock Acceleration 

The acceleration of thermal particles, or the (re)acceleration 
of ambient Galactic CRs, to relativistic energies through DSA 
can be described by the following diffusion-advection equa- 
tion written in the shock rest-frame. 



df 

[u{x)-VAix)] — 



d_ 
dx 



D{x,p) 



dl 
dx 



- Q{x,p) 

d[u{x) -VA{x)] pdf 
dx 3 dp 



(10) 



The lab-frame density precursor is present in Figure [T]but it is much less 
pronounced than the others. The reason for this can be seen by compaiing 
p{x) and V{x). The density is given by p{x)/po = l/U{x) and when the 
shock is modified, U(x) < 1 and the incoming fluid is slowed by a fraction 
1 — U{x). When the constant wind is taken into account, the lab-frame 
plasma speed is V{x) = Vys(x) — U(x)(Vys ~ Kvind)- In this particular 
example using parameters that fit SNR J1713, U(x) ranges from 0.7 to 1.0 at 
most, accounting for the small variation in p(x). The plasma speed, however, 
is y(x)/K,i„d = [1 - U{x)\{VYs/V^ind) + U{x) and Fps/Kvind - 
200 — 300, making the speed precursor much more prominent. 

' The shaip density and temperature features at the CD are artifacts of the 
initial conditions for the hydro simulation (see the feature where the ejecta 
joins the ISM material in the dashed red curve in the top panel of Figure [T] 
and the assumption of sph erical symmetry . Similar features are present in 
self-similar solutions (e.g., IChevaliei1 [19821) . The amount of material in the 
sharp density region is generally small compared to the total in the shocked 
circumstellar medium and does not bias the total thermal emission in any 
noticeable way. 




4 6 8 

radius [pc] 

Fig. 1 . — Lab-frame snapshots of the density (p), total pressure (Ptot = 
^CR + Pg)^ g3s to total pressure ratio (Pg / Ptot), temperature (T), flow 
speed (V), and magnetic field (B) as a function of radius from the explosion 
center at equally spaced times after the explosion. The parameters used are 
those of Model II in Table[T] The red dotted lines show the initial conditions 
of the model at to = 20 yr. The parameter for the FEB is /fob =0.2 and the 
precursor is clearly evident in all of the profiles except for the density. The 
density precursor is present but is much less pronounced as explained in the 
text. The CD shows up as a sharp density drop in the top panel and the RS is 
present but no attempt has been done to accurately resolve it. The behavior 
of Pg/ Ptot comes about because the thermal pressure in the pre-SN wind is 
very low due to the low temperature and small gas density which drops with 
radius, such that Ptot is always dominated by the CR pressure Pcr in the 
precursor. At the subshock, Pg/ Ptot jumps from ~ 0.003 to ~ 0.7 mainly 
due to the strong shock heating of the upstream gas to T > 10® K, and 
becomes unity behind the CD since particle acceleration is not considered on 
the reverse shock side in this work. 



where the diffusion coefficient D{x,p) depends both on the 
momentum of the particles and the position in the precur- 
sor, and Q{x,p) is the injection rate from the shock-heated 
plasma. The intrinsic assumptions that the phase-space distri- 
bution of the accelerated particles f{x,p) is nearly isotropic 
in the fluid frame, and that the shock is non-relativistic are 
used in this equation. The diffusion coefficient should vary 
with x as the local i?-field changes due to MFA in the CR 
precursor. As an option, we can parameterize D{x,p) in the 



following form: 



D{x,p)=Do 



3fiG 
W) 



P 



(- 



Po 



(11) 



where pq = 1 GeV/c unless otherwise specified, /3 = v{p)/c 
with v{p) being the particle velocity, and Dq is a normaliza- 
tion factor. Here, B{x) is the amplified magnetic field in the 
precursor region. For Bohm diffusion, for example, we have 
a = 1 and Dq = poc^/[3e{3fiG)]. 

Alternatively, we can compute the diffusion coefficient 
by extending a quasi-linear theory of wave generation (e.g., 
iBlandford & Eichler 1987) to the high-amplitude regime as 
a first approximation. In this case, D{x,p) is calculated 
from the self-generated wave spectrum with a power-law form 
W{x^k) = WQ(x){kQ/k)^~°', where k is the wavenum- 
ber and is related to the particle momentum through gyro- 
resonance, i.e., fcies = eSB{x)/{pc). Note that this res- 
onance relation is space dependent due to MFA in the CR 
precursor The normalization factor Wo{x) can be obtained 
self-consistently using the pressure from magnetic turbu- 
lence Pju through equation ^ and the relation Piu{x) = 
^ / W{x, k)dk for Alfven waves. The resonant scattering 
mean-free-path is given by Arcs = {^ / ''^)Pw{x) / [k'^W {x , k)]. 
Given these assumptions, we can readily obtain the diffusion 
coefficient as Arcsw/3, 



D{x,p) 



37Te5B(x) 
2p, 



\Z!) ' 



3ir{l-a}e5B{x} \p, 



ifa = 1 
if a < 1 



It can be shown that the distribution of accelerated particles 
at the subsh ock can be well approximated by the following so- 
lution (e.g., lBlasill200l fAmato & Blasll2005l: Ic^ioli et al.l 
l20T0ah : 



frip) - 



35t, 



g— CXP 



47rpj 



iiij 



StotU{p) - 1 

P dp' 3StotU{p') 
p,^^ p' StotUip') - 1 



' Pmin " 



dp' 3StotU{p') 



P' StotUip') - 1 



(13) 



Here we have adopted th e therm al injection model de- 
scribed in e.g. iBlasi et al] (120051) , such that Q{x,p) = 
[r]niui / {A'Kpf^^)]S{x)S{p — pinj) is the injection rate, ?y = 



3v^('^«- 



i)xr 



inj" 



-Xi„ 



is the fraction of thermal particles 
being injected into DSA, and xinj is defined such that the 
injection momentum in the shock frame is pinj = (xinj — 
U2/c)^2mpkT2 (e.g.. ICaprioh et al.ll2010bh . The 'cr' su- 
perscript on ,f"{p) indicates that only superthermal particles 
are included in equation (fTsT l. The first term in the braces in 
equation ( fT3] l describes the acceleration of the injected par- 
ticles from the downstream thermal pool, while the second 
term describes the acceleration of the pre-existing CRs in the 
ISM, with an ambient spectrum fooip)- Here, for simplicity, 
we have assumed that the minimum momentum Pmin = Pinj 
for the pre-existing CRs. To produce a full spectrum of all 
shocked particles, /i(p), we add a thermal component to 
/" (p) with a density and temperature that are consistent with 



the NL DSA calculation, producing a spectrum similar to the 
red solid curves shown in Figure |3] 

The 'effective' compression ratios S'tot and Ssuh experi- 
enced by the streaming particles scattering off the magnetic 
turbulence are smaller than those felt by the background fluid 
if the scattering centers are moving away from the shock up- 
stream in the plasma frame. We assume that the magnetic 
scattering centers have a local speed, relative to the plasma, 
of 



va{x)-^ 




- 1 



/alf 



/. 



alf 



(14) 



in the direction parallel to the shock normal and away from the 
shock. The model parameter /aif G [0, 1], allows us to vary 
va{x) linearly between Bq/ y/4TTp{x) and B{x)/ ^/A'Kp{x). 
When /alf = 0, v^ corresponds exactly to the definition of the 
Alfven speed, which is believed to represent well the veloc- 
ity of the magnetic turbulence in the limit of small-amplitude 
5B <C i?o- As the self-generated turbulence deviates from the 
small-amplitude regime, however, the wave speed may also 
deviate from the Alfven speed. We use the linear parametriza- 
tion in equation (fl4] | as a simple model to reflect this transi- 
tion. Similar approac hes have been used previously by sev- 
(e.g. iBell & LuceklbOOli IVladimirov et al.l 



C\j\ eral other groups (e.g. 



l2008Hzir^shvili & Ptuskinll2008HCaprioli et al.ll2008l) . We 
note that this recipe can be modified as theories and/or sim- 
ulations will improve. With this definition of va{x), we can 
obtain: 



Up - VAfl 
U2 + VA,2 

■ Rtot 



1 - VAn 



Ss- 



;ub ' 



1 + Va,o[^ + {B2IB0 - l)./alf]/R^ 
U2 + VA,2 

1 - Va.o[^ + (Si/Bo - l)/aif](7r'^'' 



R 



■sub 



1 + VaA^ + {B2/B0 - l)/alf]v^ 



(15) 



where Va = va/uq is the dimensionless wave speed. The 
general effect of using the amplified magnetic field instead of 
the background field in equation (fT4l i is a larger wave velocity 
and a decrease in the effective compression ratios, and hence 
a less modified shock structure and reduced acceleration effi- 
ciency. This effect is especially notable for cases when MFA 
is efficient. The quantity U{p) in equation (fT3T l roughly indi- 
cates the flow speed experienced by a particle with momen- 
tum p at one diffusion length upstream from the subshock, 
and is given by: 



U{p)^ 



dx 



[U{x)^Va{x)] 

Xp{x,p) 



exp 



dx 



,Uix')-VAix') 



cp(a;',p) 



(16) 



with Xp{x,p) = 3SsnhD{x,p)/[{Ssuh - l)qip)uo] and 
g(p) = -ain[/r(p)]/51n(p). 



To determine the field compression at t he subshock, we fol- 
low previous work (e.g., lRevnoldslll998h and assume that the 
_B-field is totally random in orientation. The downstream field 

is then given by B2 ~ Bi -^/l/S + 2/3i?g^,^, where Bi is the 
field strength immediately upstream of the subshock. 

When the above elements are combined, we obtain the su- 
perthermal particle distribution in the precursor at any posi- 
tion x: 



rix,p) = /rb)cxp 



dx' 



,Ui x')-VA{x') 

Xp{x',p) 



(17) 



To obtain the full spectrum, f{x,p), a thermal portion has 
been added at each x with a density and temperature consis- 
tent with the precursor values at that position in a similar way 
to f"{p) at the subshock. 

2.2.1. Maximum particle momenta 

The determination of the maximum momenta of the ac- 
celerated protons and electrons is done self-consistently us- 
ing the precursor information and diffusion coefficient. For 
protons, Pmax is obtained by either equating the acceleration 
time-scale with the remnant age (age-limited), or by equat- 
ing the upstream diffusion length with a 'free escape bound- 
ary' (escape-limited) set in front of the shock at a distance 
Lfob = /fcbi?SNR(i), where /fob is typically set between 0.1 
and 0.2. Whichever limit gives the lowest Pmax is used. The 
acceleration time-scale is approximated as (e.g.. lDrurvlll983l : 
iBerezhko & Ellisoniil999i) : 



U 



Uq - U2 Jp 



dp fD„ip) , D^ip) 



P 



Uq 



U2 



(18) 



To obtain the escape-limited Pmax, we equate a diffusion 
length averaged over the precursor, i.e., 

{D{x,Pnje,x)/u{x)) ^ / [D{x, Pmi,x) /u{x)]dx / Lf^h , 

(19) 
with the free escape boundary distance. Note that in the case 
when the diffusion coefficient is determined from the self- 
generated wave spectrum, D{x,p) is dependent on Pmax as 
shown in equation ( fT2b . Therefore, Pmax has to be obtained 
through recursive steps in the code. 

For electrons, Pmax can also be limited by radiation losses, 
mainly from synchrotron and inverse-Compton losses. The 
loss tim e scale of a relativistic e lectron, t\oss — ill, is given 
by (i.e.. lSchUckeiser et al.ll2010l) : 



AutcUbI 



Wd 



WriL 



It UB{i^ + ^h) 



(20) 



where we define Ub ~ ^2/8"'' 7 is the Lorentz factor of 
the electrons, A^phot is the number of seed photon fields, and 
Ik.i = Q.b'ira^c^ /UbT and Wi are the critical Klein-Nishina 
Lorentz factor and energy density of the i-th seed photon field, 
respectively. For simplicity, only the CMB photon field is 
considered for the models in this paper We then obtain the 
loss-limited Pmax by equating the loss time-scale with the ac- 
celeration time-scale. When radiation loss is slow relative to 
acceleration, the Pmax of electrons is the same as for protons. 



2.2.2. Escaping cosmic rays 

The final process that needs to be accounted for to obtain 
a self-consistent solution for the CR production is the energy 
flux of particles that escape upstream from the shock at the 
FEB. When Pmax is set by escape, the highest energy CRs 
leave the FEB with a flux determined by energy conservation 
and is given by: 



F,. 



1 3 , 7ff n , 7CR,0 n 

^ 7s - 1 ■ 7CR,o - 1 



1 3 , 7s n , 7CR,2 
0P2W2 + —rPg.2U2 + 

2 7s - 1 7CR,2 - 1 



-FCR,2W2 + Fy, 



(21) 

In equation (I2TI 1. the magnetic energy flux downstream from 
the shock, in the limit va,2 ^ W2, can be approximated 

as F„,2 ~ 'Su2P^,2 = 3u2P^,i(l/3 + ^/'^RlJ- The 
shape of the escaping CR distribution is determined as in 
lEllison&Bvkovl(l201lb . 

Once the above equations are formulated, the NL DSA 
problem is solved by making an initial guess for Ui and then 
using this to calculate the quantities U{x), ^crI^:) and hence 
U{p), B{x), D{x,p), Pmax and fi{p). The self-consistent 
result, where all of the above quantities are mutually con- 
sistent and the system of equations is internally satisfied, is 
determined through i teration using the scheme described in 
lCapriolietaTI(l2010al) . 

2.2.3. NEI and temperature equilibration 

Typical SNRs have densities that are low enough that the 
ionizaton state of the shocked gas can differ significantly from 
equilibrium. The approach to equilibrium depends on the 
electron and ion densities [Ue, rii), the electron temperature 
(Te), and the time for which the shocked gas has been ion- 
izing. Since these quantities are also affected by the dif - 
fusive shock acceleration process (e.g.. lEllison et al.l 120071) . 
we have incorporated a calculation into the model that fol- 
lows the nonequilibrium ionizati on of the shocked gas (i.e., 
iPatnaude. Ellison. & Sland 120091) . and couples the resultant 
ionization vectors, densities, and electron temperatures to a 
plasma emissivity code (i.e., Patnaudeet al. 2010). The ion- 
ization structure of the shock-heated gas at a particular dis- 
tance behind the shock in a SNR is determined by the electron 
density rig, the electron temperature Tg, and the ionization and 
recombination rates for each ion of interest. This structure is 
calculated by solving the collisional ionization equations in 
a Lagra ngian gas element beh ind the shock as described in 
detail in lPatnaude et alj (120091) . 

Since the hydro simulation does not specify the individual 
electron and proton temperatures, a temperature equilibration 
model for the postshock region must be assumed. Two ex- 
treme cases for electron heating are clear: instant equilibra- 
tion with the proton temperature, and eq uilibration thr ough 
Coulomb collisions (i.e., equation 5-30 in lSpitzeJll962l) . For 
both cases we include the effects of adiabatic expansion on the 
temperature components. While all of the examples shown in 
this paper assume that the electrons are heated by Coulomb 
collisions, the cr-hydro-NEI code allows for a continuous 
range of heating modes between the two extremes. Examples 
wit h faster heating modes are shown in lEllison et al.l (12007b 
and lPatnaude et all ( 120091) . 

2.2.4. Ejfect of ion-neutral damping in a partially -ionized medium 



It is known that Alfven wave evanescence due to ion- 
neutral collisions in a partially-ionized precursor can lead 
to wave damping and hence a reduction of trapping power 
for the highest energy particles. Recent gamma-ray obser- 
vations by Fermi, combined with previous TeV observations 
by ground-based Cherenkov telescopes, have revealed high- 
luminosity, broken power-law gamma-ray spectra from nearly 
a dozen Galactic SNRs which are known to be interacting di- 
rectly with adjacent molecular clouds (e.g. lAbdo et al.ll2009l 
l2010allbh . It has been postulated immediately after the dis- 
covery that ion-neutral damping can naturally give rise to a 
spectral break in the accelerated particles for fast molecular 
shocks, which may prov ide an explanation for these peculiar 
gamm a-ray spectra (e.g.. lMaIkov et al.ll201 ItlUchivama et alj 
l2OT0h . 

In order to model shocks running into a partially-ionized 
material, we have included phenomenologica l features fol- 
lowing previous stud i es by iDrurv et all (119961) : iBvkov et al.l 
(l2000l) : lMalkov et alj(l201 lb and references therein. The pres- 
ence of a spectral break in momentum, pbr, originates from 
the fact that particles with p > pbr are only weakly scat- 
tered by non-resonant small-scale waves. The diffusion co- 
efficient is hence modified in the precursor region due to a 
transition of the wave spectrum below and above pbi ■ When 
ion-neutral damping is exp ected to be i mport ant, we fol- 
low the parameterization in iBvkov et al.l (120001) and impose 
a momentum break for the diffusion coefficient at pbr such 
that D oc {p/phr)" for p < pbr and D ex {p/phiY for 
P > Pbr, respectively. The momentum pbr can be estimated 
by equating the gyro-radius to the reciprocal of the evanes- 
cence wavenumber Fin = 2TTVa/i'in, where z/jn is t h e ion- 
neutral collision frequency. Following iMaUcov et al.l (1201 lb 
and references therein, we write: 



Pbr 



10{T, 



-0.4„ 



'Bl) rupc , 



(22) 



where the upstream quantities, no,i, are the number densities 
of neutral and ionized particles respectively, T4 = T/IC* K, 
and i?p = B/1 jiG. Assuming that the pitch-angle distribu- 
tion of the particles below the break momentum remains lo- 
cally isotropic, and that the underlying spectrum extends to 
arbitrarily high energy, the spectrum steepens by one power- 
law index atpbr, i-e., f{p) -^ f{p){Phr/p) forp > pbr- 

In addition, wave-damping due to ion-neutral collisions can 
result in a high-energy cut-off by imposing limits on the max- 
imum energy ac hievable by the accelerated p articles through 
DS A. We follow lDrurv. Duffv. & Ki^ (Il996h . who estimates 
this maximum energy to be 

, = ul.T^'^-^n-'n'l-'^CR.-i GeV , (23) 



E,, 



where uqj = uo/lO'^km s~^, and ^cR.-i ~ ^cr/0.1 evalu- 
ated at the subshock. This puts an extra constraint on Pmax in 
addition to the factors described in the previous section, when 
the upstream medium is only partially ionized. We reserve a 
detailed discussion of ion-neutral damping for future work. 

2.2.5. Cases with a pre-SN wind 

In the above equations it was implicitly assumed that the 
shock was propagating into a uniform medium as might be 
the case for thermonuclear. Type la supemovae. For core- 
collapse SNe, it is expected that pre-SN winds will produce 
spatial gradients for the pre-shock quantities like the back- 
ground _B-field, gas density, and gas pressure, which are de- 
pendent on the properties of the wind such as the mass loss 
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Fig. 2. — Evolution of a SNR with test-particle DSA (xinj = 4.7; solid, 
black curves. Model III) and efficient DSA [Xinj = 3.6; dashed (red) curves 
Model IV. and dot-dashed (blue) curves Model V]. Panel A shows the FS 
radius, panel B shows the FS speed, panel C shows the FS total compression 
ratio, panel D shows the inagnetic field, B2, immediately downstream (DS) 
from the FS (the curve marked Bq is the unshocked field for all three models), 
panel E shows the DSA efficiency at any given time, and panel F shows the 
dels all h 

0.1 cm~'^. Bo = 3 fiG (shown as a dotted line in 
panel D), To = 10* K, isNR = 1000 yr, and /^if = 0. The starting time 
for the simulation is ip = 30 yr and the transition froin age-liinited to size- 
hmited acceleration occurs smoothly with no noticeable change in the SNR 
evolution at t < 40 yr 



rate of the progenitor star. These gradients have to be taken 
into account when the recipe described in the previous sec- 
tions is applied, especially when we calculate the profiles 
of various quantities in the CR precursor. In these cases, 
we modify, for example, equations (|2]) and (0]) by replacing 
the constant background values Pgo and Bq with the space- 
dependent values in the pre-shock wind, and similarly for the 
gas density and velocity. Throughout this paper, we assume 
an isothermal and constant velocity wind whenever one is ap- 
plied. 

3. RESULTS 

We first show in Section [TT] how efficient DSA can mod- 
ify the evolution of a SNR with a simple example of a SNR 



8 



exploding in a homogeneous environment typical of Type la 
supemovae. Then, since the generalizations described above 
closely follow the work of Blasi and co-workers, we com- 
pare our results to recent work from that group both with and 
without MFA in Section [J!2l It is important to note that al- 
ternative descriptions of NL DSA exist and in Section 13.21 
we compare our resu l ts to t he Monte Carlo ca l culation s of 
lEUison & Vladimirovl (120081) and lCaprioU et al.l (l2010bh . In 
Section [331 we give a detailed discussion of SNR RX J1713.7- 
3946 and compare our results against those of lEUison et al.l 
(|2012). 

3.1. Modified evolution ofSNRs with self-generated 
magnetic turbulence 

In this section we compare the time evolution of a SNR un- 
dergoing efficient DSA with MFA to a test-particle (TP) case. 
The solid (black) curves in Figure |2] are for the TP example 
(Model III, Xinj =4.7) and the dashed (red) and dot-dashed 
(blue) curves have efficient acceleration (xinj = 3.6). For the 
dashed (red) curves, the wave damping factor /damp = 0.1 
(Model IV), while for the dot-dashed (blue) curves /damp = 
0.99, suppressing MFA (Model V). As expected, efficient 
DSA results in a slower FS speed and a smaller SNR radius at 
a given age than for the TP case. The proton Pmax is increased 
dramatically for the NL case and is > 100 times the TP Pmax 
for /damp = 0.1. When /damp approaches 1 and MFA is sup- 
pressed, the changes between the TP and NL results are less 
pronounced but the proton Pmax is still ^ 10 times greater 
in the NL case. As shown in Tabled we note that the Pmax 
of electrons is limited by fast synchrotron cooling for Model 
IV due to the highly amplified downstream B-field B2 from 
efficient MFA. On the other hand. Model III and V produce 
accelerated electrons with Pmax's essentially the same as the 
protons. In these cases, MFA is inefficient due to a low DSA 
efficiency (TP) and a high wave damping rate respectively, 
and the resulted B2 is not large enough to have any impor- 
tance on determining the electron cutoff's. Also, in this ex- 
ample, since Sq = 3 jiG is relatively small, the magnetic field 
has little influence on the SNR evolution as indicated by the 
fact that the dashed and dot-dashed curves for _Rfs and Vps 
are almost indistinguishable. 

3.2. Comparison with \Caprioli et al\^010cm ) 

In the top panel of Figure |3] w e show a proton spectrum 
presented in ICaprioli et al.l (l2010al) (dashed curve) compared 
against our static comparison model A (solid curve). The pa- 
rameters for the cr-hydro-NEI model are shown in Table |3] 
and both of these models include MFA. All spectra shown in 
Figure |3] are static examples with no SNR evolution included 
and fi{p) is measured at the subshock. The superthermal 
portion of /i (p) below the high-energy turnover is extremely 
well matched in both shape and absolute normalization, as is 
the value of Pmax- We note that we determine the shape o f 
the turnover at Pmax differentl y thanjCa prioH et al .l (12010a 
Follo v ying previous work (e.g., lElhson e t al. 2004t lLee et al.l 
l2008t lEUison et al.l 120 12h . we parameterize the turnover by 
multiplying f{x,p) by an exponential factor, i.e.. 
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f{x,p) -^ f{x,p)exp 



1 



P 



(24) 



where acut is a free parameter and /3cut = 1 or acut depend- 
ing on the particular application. Results for different values 
of ttcut are shown and labeled in the top panel. In all of the 
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Fig. 3. — Top panel - Comparison between our static Model A (liglit-weiglit 
colored curves) and tlie efficient acceleration example in ICaprioli et al.l 
)2010al) (heavy black dashed curve). Both of these models include MFA. 
Curves of different colors con'espond to different cicut assumed: Ocut = 1-5 
(red solid), ctcut = 0.75 (blue dashed) and ctcut = 3.0 (green dash-dotted). 
Middle panel - Comparison between our static Model B (solid red curve) and 
the benchmark model in Caprioli et al. (2010b) (dashed curve - semi-analytic 
model). The dotted curve in the middle panel is the Monte Carlo model from 
jCaDiioli et al. 2010^). These models do not include MFA. Bottom panel 
- Comparison between cr-hydro-NEI (model C) (solid red curve) and the 
Monte Carlo result from|Ellison & Vladimirov (2008). Both of these models 
include MFA. The v ertical solid lines indicate the pmax obtained by our static 
models (see Section [2.2. It In the top panel, results with different cut-off in- 
dices (acut in equation l24l are shown as indicated. In the middle and bottom 
panels, the cr-hydro-NEI model assumes Xinj =3.1 and 3.5 respectively, 
and Qcut = 1.5 for both. In all cr-hydro-NEI lesuits, /3cut = «cut- 



cases shown in Figure[3] /3cut = ctcut- The vertical solid lines 
indicate the Pmax obtained from the cr-hydro-NEI model. 

While we coul d have adapted the procedure described 
in 'Caprioli et al.l (l2010ah (or other models such as 
jZirakashvili & PtuskinI 120081) to determine the shape of 
the turnover, we feel the fundamental uncertainty in this 
shape is currently bes t described by free parameters (see 
lEUison & Bvkovir201 IL for a fuller discussion). The CRs that 
escape upstream from the shock will have a highly anisotropic 
distribution and are certain to generate turbulence that will 
determine the nature of the F EB and the shape of the escaping 
population (see, for example. iBykov et al. 201 lb). Since the 
shapes of the trapped and escaping CR distributions, at the 
highest accelerated energies, are critical for modeling both 
X-ray synchrotron emission and GeV-TeV 7-ray emission, 
we feel it is important to have a parameterized model that can 
be compared to observations until a more definitive theory of 



self-generated, long-wavelength turbulence is developed. 

A clear difference in the results displayed in the top panel 
of Figure [3] is that we include the thermal population in 
flip)- The thermal part of fi{p) is calculated assuming it is 
Maxwell-Boltzmann with the shocked temperature and den- 
sity that are self-consistently determined by the DSA calcula- 
tion. A similar technique is us ed for the semi-analytic calcula- 
tion in ICaprioli et al.l (l2010bh and the full spectrum is sh own 
in the middle panel for both the (ICaprioli et al.l 1201 Obi) and 
cr-liydm-NEI models^ The models in the top panel include 
MFA while those in the middle panel do not. 

The examples shown in the top two panels of Figure |3] and 
other tests we have performed but do not show, convince us 
that we have successfully included the static NL DSA calcu- 
lation developed by P. Blasi and co-workers in the cr-hydro- 
NEI code. Ou r implem entation is not identical to that used in 
ICaprioli et al.l (l2010al) because we have chosen to treat some 
processes, such as CR escape and the calculation of the dif- 
fusion coefficient, in a different fashion. However, a suitable 
choice of parameters can yield essentially identical results as 
Figure [^demonstrates. 

This in not the case for the Monte Carlo comparisons shown 
in the middle and bottom panels of Figure|3] The Monte Carlo 
simulation makes fundamentally different assumptions from 
the semi-analytic calculation, particularly for how injection 
is treated and how the highest energy CRs escape from the 
shock. These differences show up in the spectra shown in 
Figure [3] In contrast to the semi-analytic model, the Monte 
Carlo simulation does not make a diffusion approximation, 
i.e., it is never assumed that the particle distribution is nearly 
isotropic in some frame. Since thermal particles are never 
nearly isotropic in the shock frame, the injection of thermal 
particles cannot be directly modeled with a diffusion equation 
such as ( [Tol l. Escaping CRs will also have a highly anisotropic 
distribution making them hard to model within a diffusion ap- 
proximation. 

For the pa r ameters used for the "benchmark case" in 
ICaprioli etalj ( 1201 Obi) , the Monte Carlo model yields a high 
injection efficiency and Xinj ~ 3.1 was chosen to match it for 
both of the semi-analytic models in the middle panel. While 
this Xiiij yields an acceleration efficiency and total compres- 
sion ratio that is very similar for the semi-analytic and Monte 
Carlo results, the spectrum in the transition region between 
thermal and relativistic CRs is quite different. Of course, the 
Monte Carlo results depend on the assumptions of the model 
which do not include a description of the complex plasma 
processes that are certain to be important in the subshock 
layer. The shapes of the spectra near Pmax are very simi- 
lar in the three cases shown in the middle panel where we 
have used acut = /?c ut = 1.5 for the cr-hydm-NEI case (see 
ICaprioli et al.l 1201 Obi for a thorough discussion of the shape 
of the turnover near Pmax)- 

In the bottom panel of Figure [3] we compare the cr- 
hydm-NEI model to t he M onte Carlo results presented in 
Ellison & Vladimirovl (120081) . An injection parameter Xinj = 
3.5 is chosen to match the total compression ratio obtained 
by the Monte Carlo simulation (i?tot ~ 9). A similar differ- 
ence in the transition region between thermal and superther- 
mal shows up in the two models as expected. Here, however. 



We note that while the semi-analytic solution determines the pressure 
accurately, the "thermal" part of f{x, p) does not have to be strictly Maxwell- 
Boltzmann in shape immediately behind the shock. This will depend on the 
equilibration time scale. 
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Fig. 4. — Broadband spectrum and the "best-fit" model for SNR RX 
J1713. 7-3946 listed in Table [T] as Model I. The different emission pro- 
cesses are: synchrotron (solid blue curve), IC (dot-dashed purple curve), 
pion-decay from trapped CRs (dashed red curve), bremsstrahlung (dotted 
black curve), and thermal X-rays (solid black curve). The dash ed black 
ojrve is the summed emission. The data is from 'Acero et al.l I I2009I) (radio), 
Tanakaetal. (2008) (Suzaku X-rays), Abdo et al. (2011;) (Fermi-LAT), and 
lAharonian et al.. (.20111) (HESS). Note that the two lowest energy Fermi-LAT 
points are upper limits. A column density of uh = 7.9 X 10^^ cm^'^ has 
been used for absorption and the inodel spectra have been multiplied by a fac- 
tor /norm = 0.65 to match the overall normalization of the observations. As 
noted in the text, the low-energy thermal X-ray emission is fully consistent 
with the Suzaku data once the instrument response is considered. 

there are clear differences in the shapes of the spectra over a 
larger fraction of the momentum range. Particularly important 
is the difference in the location of the minimum in the p^f{p) 
plot; the Monte Carlo minimum is at a momentum ~ 5 times 
as high as the semi-analytic result. 

3.3. Pre-SN wind model of SNR RX J1713. 7-3946 revisited 

We now examine how the modeling of SNR J1713is modi- 
fied using the generalized cr-hydro-NEI code compared with 
a previous p r e-SN wind model using the code described in 
lEllison et al.l (120121) and references therein. A 'best fit' result 
using our generalized code is shown in figures |4]and|5](Model 
I in Table [Tji and it is clear, w hen Figure |4]is compared to fig- 
ure 2 in lEllison et al.l (l20I2b . that our generaliz ed model can 
yield a fit comparable in quality to that shown in Ellison et 
(20121. The most important conclusion from lEllison et al.l 
( I2OI2I) . that inverse-Compton emission from electrons scat- 
tering off the CMB photo ns dominates the G eV-TeV emis- 
sion, remains robust. As in lEllison et al.l (120 1 2b . the emission 
shown in Figured is the integrated emission from the inter- 
action region between the CD and the FS and the shocked 
electron temperature is determined assuming the electrons are 
initially cold [i.e., T^ = {me/mp)Tp\ and then heated by 
Coulomb collisions, as discussed in Section l2.2.3l In fact, in 
this case, the dominance of IC over pion-decay at GeV-TeV 
energies doesn't depend on how electrons are heated although 
specific X-ray line ratios do. The fitting paramet ers used here 
for Fig ure |4] are close enough to those used in lEllison et al.l 
( I2OIOI) so the detailed discussion of the electron equilibration 
in that paper appli es here. We refer the reader to text around 
Figures 1 and 2 in lEllison et al.l (120101) for a full discussion. 
The emission from the region between the reverse shock and 
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Fig. 5.— Evolution ofModelIforSNRJ1713. Panels A, B, and C show the 
FS radius, speed, and Mach numbers, respectively. In panel D the solid curve 
is the total compression ratio, Rtot , and the dashed curve is Ksub ■ The mag- 
netic fields corresponding to far upstream. Bo , immediately upstream, Bi . 
and immediately downstream from the FS, B2. are shown in panel E. Panel F 
shows the instantaneous DS A efficiency and panel G shows the fraction, /sn , 
of supernova explosion energy, -EsN. put into all CRs (black curve) and only 
escaping CRs (red curve). The maximum particle momentum shown in panel 
H is essentially the same for protons and electrons since radiation losses are 
small for this example. The transition from an age-limited to a size-limited 
Pmax occurs at t < 30 yr. 



the CD is ignored and the pion-decay emission from escap- 
ing CRs is not shown in Figure |4] because we believe this 
w ill be insignificant f or SNR J 17 13 for the reasons detailed 

in lEllisonetaI1(l2012h. 

We note that, as in lEllison et al.l ( l2010ll20T2l) . we have cho- 
sen parameters for SNR J1713 that produce some thermal X- 
ray emission above the observations at the low-energy end 
of the Suzaku energy range. While parameters could have 
been used that would produce no thermal emission above the 
Suzaku observation, it must be kept in mind that the low- 
energy range is uncertain due to interstellar absorption and 
the Suzaku instrument response. The model we show in Fig- 
ure m is fully con sistent with the Suzaku observations (see 
et al.ll2010l for a full discussion). Any modification 



son I 



IeIH 

of the parameters to lower the low-energy thermal emission 
would strengthen the conclusion that IC dominates the GeV- 
Te V emission. This is also the case if any thermal X-ray emis- 
sion from the RS is included. 



It is important to note that, because of the number of pa- 
rameters in the cr-hydro-NEI model, there is a range in each 
parameter (typically ±20%) that allows a fit that is essentially 
equivalent to that shown in figure|4]when other parameters are 
adjusted accordingly. To illustrate this, we list the parameters 
for Model II in Table[T]that produce a plot (not shown) almost 
indistinguishable from Model I. Despite this freedom, there is 
no set of parameters that allows pion-decay to dominate IC at 
GeV-TeV energies as long as the constraint from the thermal 
X-rays emission lines is considered. The low thermal X-ray 
luminosity forces the ambient density to be low and this re- 
sults in IC dominance at 7-ray energies. While t he fit to SNR 
J1713 shown here is similar to that shown in lEUison et al.l 
( I2OI2I) . there are differences in the fitting parameters stem- 
ming from the way MFA, the CR diffusion coefficient, Pmax, 
the shock precursor, and other aspects of NL DSA are calcu- 
lated. 

In the model shown in figure 2 of Ellison et al.l (120 1 21) . the 
acceleration efficiency was set at ^dsa = 25% and the in- 
jection parameter, Xi^j, was determined as the SNR evolved 
to produce this percentage^ Here, we fix xinj = 3.75 and 
the acceleration efficiency is determined directly from this PI 
In both cases, either ^dsa or Xinj was chosen to obtain a 
good fit to the observations. With xinj = 3.75, the in- 
stantaneous acceleration efficiency varies smoothly between 
£dsa — 25% early in the evolution (i.e., tsNR — 200 yr), and 
£dsa — 35% at the end of the run (isNR = 1630 yr). At 
isNR = 1630 yr, the fraction of SN explosion energy, Ea,-^, 
that has been placed into CRs, bo th trapped and escaping, is 
- 16%. In the lEllison"eran (|2072|) fit. ^ 15% of £;sn wasput 
into CRs by 1630 yr We emphasize that although IC emis- 
sion from electrons dominates the broadband emission, the 
SNR forward shock is accelerating CR ions with the efficien- 
cies just quoted and the CR electrons carry a much smaller 
frac ti on of the ene r gy. 

In lElUson et alj ( 120 1 21) an 'ad hoc' amplification factor, 
Samp = 8.5, was used where the compressed magnetic field 
at the shock was simply multiplied by a factor chosen to ob- 
tain a fit to the observations. Here, the MFA is calculated 
in a more self-consistent manner as described in Section |2] 
In both cases, however, parameters are chosen to produce a 
shocked magnetic field of ^ 10 ^G, as required to fit the ob- 
servations. The shocked magnetic field depends on a num- 
ber of param eters including cTwind = B{rY /['^'Kp{r)V^-^^^ 
dChevalier & Luo 1994i) and the amount of damping parame- 
terized with /damp- For Model I, CTwind = 0.004 and /damp = 

0.1, while for Model II, CT„ind = 0.008 and /damp = 0.3. 
As shown in panel E of Figure |5] the pre-SN wind magnetic 
field, So is increase d by ^ 100 fold at tgNR = 1630 yr As 
mentioned in Ellison et al.l ( 120121) . the low magnetic field val- 
ues we find are quite different from the high-field estimates 
(~ ImG) obtained for SNR J 17 13 from r apid time varia- 
tions and sharp X-ray edges in filaments (e.g.. lUchiyama et aP 
l2007h . However, other estim ates yield lo wer values (seerefer- 
ences in lElUson & Vladimir ov 2008), and lBykov etaH (120081) 
provide an alternative explanation for rapid time variations 
that does not require such large fields. Furthermore, our esti- 



' The acceleration efficiency, fDSA. is the percentage of far upstream en- 
ergy flux crossing the shock at any given time (in the shock rest frame) that 
ends up in superthermal particles. 

'" This change is implemented mainly to speed computation. In our current 
cr-hydro-NEI code, runs with fixed ^dsa take significantly longer to run than 
those with fixed Xinj ■ 
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mate is for the integrated magnetic field strength in the CD-FS 
interaction region and local values in small regions could be 
much higher than our average value. 

An important parameter for DSA in general, and a crit- 
ical one for determining the relative importance of IC vs. 
pion-decay at GeV-TeV energies, is the electron to proton ra- 
tio at relativistic energies, A'op. This ratio has not yet been 
reliably determined from theory or plasma simulations and 
a wide range of values, i.e., 0.02 > Kc^> 10~^, have been 
_used to fit SN R J1713 (s ee, for example, Katz & Waxman) 

20101 



2008: Morlino et all 12009 : Zirakashvili & Aharonian 



Ellison et al.ll2012l and references therein). Values of K 



op 



0.01 allow pion-decay to dominate the gamma-ray emission 
but are not consistent with the constraint imposed by the weak 
thermal X-ray luminosity. Here, as in our previous work, we 
find that A'cp ~ 0.01 provides an excellent fit to the broad- 
band emission and that values A'cp ^ 0.01 can be excluded 



The semi-analytic DSA model used in lElhson et al.l (120 12h 
assumed the CR diffusion coefficient was a rapidly increasing 
function of momentum without explicitly specifying the mo- 
mentum or spatial dependence. Our generalized model used 
to produce the fit shown in Figure |4] uses the explicit form 
given in equation (fTTI) with a = 1 and a spatial dependence 
determined by B{x). 

The determi nation of Pma x is als o more consistently done 
here than in lEllison et al.l (120 12h . as described in Sec- 
tion 12.2.11 However, when the size of the shock determines 
Pmax, as is the case for Figure IH there still must be a free pa- 
rameter, /fcb, relating the finite shock radius to Pmax- Here, 
■ff ph = 0.2 whi l e ff^h = 0.1 was used for the fit in figure 2 
of lEllison et al.l (120121) . There are a number of other minor 
differences as seen by comparing Model A in Tables 1-4 in 
lEllison et al.1 ( 120121) to Model 1 in Tables[T]and|2]here. 

In Figure |5] we show the evolution of some parameters over 
the Ufetime of SNR J1713. At isNR = 1630yr, the magnetic 
field in the pre-SN wind has dropped to Bo ^ 10^^ G and 
the field just behind the FS, B2, is ^ 50 times larger The 
instantaneous acceleration efficiency, f dsa (panel F), and the 
FS compression ratio i?tot (panel D) increase noticeably but 
over the brief the age of this remnant, Pmax is essentially un- 
changed at 



10^ TTipC. 



4. CONCLUSIONS 

We have performed an extensive generalization of a code 
{cr-hydm-NEI) that models the evolution of a spherically 
symmetric SNR undergoing efficient, nonlinear DSA. The 
semi-analytic model for DSA that we us e is based on previous 
work by P. Blasi and co-workers (e.g., ICaprioli et al.ll2010ah 
and here we present a detailed discussion of the DSA mech- 
anism and how it is incorporated in the hydrodynamic model 
of an evolving SNR. 

Relativistic particle populations are the signature of violent, 
nonthermal activity in astrophysical sources. These sources 
are ubiquitous and many of these populations are undoubt- 
edly produced by collisionless shocks. Supernova remnants 
offer a promising platform for studying the shock acceleration 
process because they are well studied in a broad wavelength 
band, there are a number of different remnants in different 
environments, it is certain that SNR shocks produce relativis- 
tic electrons, and because it is likely that SNRs are the main 
sources of CR ions up to at least ~ 10^^ eV. 

As is clear from our discussion in Section |2] efficient DSA 
is intrinsically complicated and there are a number of param- 
eters required to properly model the NL coupUngs between 



the SNR hydrodynamics, MFA, CR acceleration, and the non- 
equilibrium ionization production of thermal X-ray emission. 
Given these complexities, our aim is to develop as complete a 
model as possible that can be used to interpret the broadband 
continuum and line emission from remnants in a reasonably 
self-consistent fashion. 

The most important generalizations included in cr-hydm- 
NEI are explicit calculations of the shock precursor structure 
and MFA. These effects are coupled to a momentum and po- 
sition dependent CR diffusion coefficient, the calculation of 
the maximum particle momentum, Pmax, and a finite Alfven 
speed for particle scatterers. Other additions include the abil- 
ity to accelerate a superthermal seed population of CRs in ad- 
dition to thermal particles and wave damping due to partially- 
ionized material. We have performed a number of tests to en- 
sure that we accuratel y reproduce the semi-analytic results of 
ICaprioli et al.l (l2010allbl) where applicable, and we have noted 
where our model of NL DSA differs from that work. We have 
also compared our results to Monte Carlo simulations and dis- 
cussed important differences of this alternate DSA model. 

There are two main features that distinguish cr-hydm-NEI 
from other SNR models. The first is that we couple the rem- 
nant dynamics to the CR acceleration and the second is that 
we self-consistently include a NEl calculation of thermal X- 
ray emission with the dynamics and CR acceleration. It turns 
out that the thermal X-ray emission is the key factor in deter- 
mining the dominant GeV-TeV emission mechanism in SNR 
J 17 13. The absence of thermal X-ray emission lines in this 
remnant forces the ambient density to be low enough that 
inverse-Compton emission from relativistic electrons domi- 
nates pion-decay from relativistic ions. Without a consis- 
tent description of the thermal X-ray emission, parameters 
can be found for SNR J 17 13 with either mechanism domi- 
nating the GeV-TeV emission (e.g.. lMorlino. Amato. & Blasil 
2009). It is important to note, however, that depending on 
the ambient conditions, it is possible that the 7-ray emission 
in other SNRs is dominated by pion-decay, as claimed for 
Tycho's SNR (i.e., iMorUno & Capriol]||201 ll) , is some near- 
equal mix of IC and pion-decay, as suggested for CTB 109 
(G10 9. 1-1.0) (Castro et al. 2012), or is even bremsstrahlung 
(e.g.. lBykovetaiJi2000l) . 

In the interest of conciseness, we have not shown examples 
here with a superthermal seed population of CRs, ones show- 
ing the effects of ion-neutral damping in partially-ionized ma- 
terial, or examples where Coulomb losses are important for 
low energy CRs. These processes are expected to be impor- 
tant for SNRs interaction with dense material and in future 
work we will use cr-hydro-NEI to model the interaction of a 
remnant FS with a dense molecular cloud. 

The authors acknowledge important discussions with A. 
Vladimirov and D. Patnaude concerning this work. D.C.E. 
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NNH04Zss001N-LTSA, 06-ATP06-21, and NNX11AE03G. 
S.N. acknowledges support from Ministry of Education, Cul- 
ture, Sports, Science and Technology (No. 23105709), Japan 
Society for the Promotion of Science (No. 19104006 and No. 
23340069), and the Global COE Program, The Next Genera- 
tion of Physics, Spun from Universality and Emergence, from 
MEXT of Japan. 
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TABLE 1 
Model Parameters - Input 





Model I " 


Model II " 


Model III 


Model IV 


Model V 


Remarks 


tSNR [yr] 


1630 


1630 


1000 


1000 


1000 


SNR age 


Xinj 


3.75 


3.75 


4.70 


3.60 


3.60 


Injection parameter 


/damp 


0.1 


0.3 


0.1 


0.1 


0.99 


Wave-damping factor 


a 


1.0 


1.0 


1.0 


1.0 


1.0 


PLindexofD(x,p) •> 


Kep 


0.01 


0.01 


0.01 


0.01 


0.01 


(e-/p) ratio at pniax " 


Qicut 


0.5 


0.5 


1.0 


1.0 


1.0 


Exp. cut-oif index "^ 


/Scut 


1 


1 


1 


1 


1 


Exp. cut-off factor "^ 


/alf 


0.25 


0.25 


0.0 


0.0 


0.0 


x-dependance of va "^ 


/fob 


0.20 


0.20 


0.20 


0.20 


0.20 


Llcb/RpS 


dAf/dt[MQyr-i] 


7.5x10-6 


7.5x10-6 


- 


- 


- 


Pre-SN Mass loss 


Kvind [kms-l] 


20 


20 


— 


— 


— 


Pre-SN wind speed 


'^wind 


0.004 


0.008 


— 


— 


— 


Wind magnetization 


-^phot 


1 


1 


— 


— 


— 


No. of seed photon fields ' 


/norm 


0.65 


0.70 


— 


— 


— 


Flux normalization ^ 



* For all models listed in this Table, MFA is included and the fraction of helium by number is taken to be /hc = 0.0977. 

^ Best-fit models for SNR 11713. These models all have Esk = 10^^ erg, A/cj = 3 Mq, tsNR = 1630 yr, dsNR = 1 kpc and T„ 

^ See e.g., equation UU . Here, a = 1 refers to Bohm diffusion, for which we set D{x, p) = pcv{p) /[3eSB(x)]. 

'^ If radiation losses are important for electrons, K^p is measured at relativistic energies but below the cutoff. 

•* See equation 1241. 

" See equation J14t . 

^phot = 1 means only the CMB photons are considered for IC emission. 
8 Values other than 1 can reflect uncertainties in the distance and/or the fraction of FS surface that accelerates CRs. 
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TABLE 2 
Model Parameters - Output t 





Model I 


Model II 


Model III 


Model IV 


Model V 


Remarks 


Ms,o 


291 


291 


251 


219 


220 


Sonic Mach number 


Ma,o 


2190 


1550 


216 


188 


189 


Alfven Mach number 


Ffs [km/s] 


4.38x10^ 


4.39x10^ 


3.78x10^ 


3.30 xl03 


3.32 xl03 


FS speed 


Rfs [pc] 


9.31 


9.30 


6.91 


6.38 


6.45 


FS radius 


Sdsa %] 


39.7 


38.8 


3.46 


83.8 


79.6 


Total DSA efficiency 


£:osc[%] 


7.60 


7.43 


3.43 


32.4 


27.9 


Efficiency for escaped CRs 


Rtot 


5.05 


5.01 


4.00 


9.20 


8.39 


Total comp ratio 


Rsub 


4.00 


3.99 


4.00 


3.63 


3.68 


Subshock comp ratio 


B2 [fJ,G] 


6.71 


6.91 


10.1 


154 


17.5 


B -field right behind FS 


Pmax (p) [rripc] 


9.35 xlO^ 


9.59 xlO^ 


1.22 xlO^ 


1.19x10^ 


1.19x10* 


Max proton momentum 


Pmax (e) [_mpC] 


9.35 xlO^ 


9.59 xlO^ 


1.22 xlO^ 


9.21 xlO^ 


1.19x10* 


Max electron momentum 


icR(oo = 0) 


0.21 


0.20 


1.17x10-* 


0.64 


0.55 


CR pressure at subshock 


(P^/Pga.)l 


0.78 


0.36 


2.60x10-2 


3.86 


4.81x10-3 


Wave/gas pres. at subshock 


XCR 

















Pressure in pre-existing CRs 



Output values at the end of the simulation (t = tsNR). 



TABLE 3 
Static Model Parameters * 



Rfs [pc] 


Ms 


Ma 


Usk [km/s] 




Bo [mG] 


no [cm 3] 


To [105 K] 


Xinj 


MFA 


Model A" 1.5 


135 


150 


5000 




5.0 


0.1 


1.0 


3.30 


yes 


Model Bl" 0.1 


30 


42 


5000 




3.0 


0.003 


20.0 


3.10 


no 


Model C-- 1.0 


30 


140 


3000 




10.0 


1.0 


7.3 


3.50 


yes 


t For all models in this Table 


, «cut = /3ci 


at = 1.5, 


/ho = 0.0, 


and 


/feb = 0.1. 


For Model A, the additional cases 


acut(= /3cut) 


= 0.75 


and 3.0 are plotted. 

" Comparison withlCaoiioh el 


;al.l(l2010al). 


















*" Compaiison with Camioli el 


tal. (2010b). 
Vladimirovi J 


















'^ Comparison with Ellison & 


20081). 
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